Diffuse wave density and directionality in anisotropic solids 
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Several general results are derived for diffuse waves in anisotropic solids, including concise expres- 
sions for the modal density per unit volume d(tj), and for the participation factor matrix G. The 
latter is a second order tensor which describes the orientational distribution of diffuse wave or re- 
verberant energy, and reduces to the identity I under isotropy. Calculations of G for a variety of 
example materials show significant deviation from I even under moderate levels of anisotropy. 
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I. INTRODUCTION 

We consider how material anisotropy effects the di- 
rectional partition of reverberant or diffuse wave energy. 
Diffuse waves in solids are the long time response when 
multiple scattering has equilibrated the energy distribu- 
tion among modes. Preferential orientation of the root 
mean square particle velocity does not arise in isotropic 
materials but is a characteristic of anisotropy. Our objec- 
tive is to describe this orientation effect and to quantify 
it in real materials. An ability to determine, directly or 
by inference, the orientational distribution of kinetic en- 
ergy density in a solid allows one to essentially "hear" 
the texture of a crystal. We will demonstrate that the 
key quantity that needs to be measured is the autocor- 
relation function, or the Green's function evaluated at 
its source. By deriving an explicit formula for the auto- 
correlation, or the admittance matrix, wc can completely 
describe the directional distribution of the diffuse wave 
energy. 

We introduce two quantities for the description of re- 
verberant energy in the presence of anisotropy: the par- 
ticipation tensor G and the modal spectral density per 
unit volume, d{uj). The two arc in fact intimately related 
as we will see. Under steady state time harmonic con- 
ditions the total energy of a body is equally divided be- 
tween potential and kinetic. The latter is / dFp|up 
where |u| is the root mean square particle displacement, 
and assuming a uniform spatial distribution, the total 
energy is E = Fpw^jup. This may be inverted to ex- 
press the mean square displacement. Let Ui = |u • e^j 
where e,;, i ~ 1,2,3 is an orthonormal triad. Since 



|u| we may write 



E 



(1) 



for i=l,2,3 (no sum) where G is a second order symmetric 
tensor satisfying 

tr G = 3. (2) 

For isotropic materials G is simply the unit matrix or 
identity (second order) tensor. Deviations from this can 



occur under three general situations: (i) If the field point 
is near a surface or boundary. This was considered in 
detail by Weaver— who found expressions for the compo- 
nents of G at a free surface in terms of simple integrals, 
see also Egle^. (ii) By analogy, G will be influenced by lo- 
cal inhomogeneity in the material, for instance if the field 
point is close to a rigid inclusion, or a void. We will not 
discuss this further here, (iii) Material anisotropy can 
also influence G. Here we consider the simplest case of a 
field point in a homogeneous material of infinite extent. 
It is expected that G displays the symmetries appropri- 
ate to the degree of anisotropy. Thus, it is characterized 
by a single parameter for materials with isotropic and 
cubic symmetries, and by two or three parameters for 
materials with lower symmetry. 

The spectral density of modes D at frequency lo in 
a volume V is D{uj) = Vd{ijj). It can be estimated as 
D = dN/dui « Vuj'^/c^ by noting the total number of 
modes scales as iV(k) « Vk^ where k = cu/c is typ- 
ical wavenumbcr. A more precise counting yields, for 
isotropic bodies, the well-known result^ 



2^ 



1 



(3) 
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where c/ and Ct are the longitudinal and transverse elastic 
wave speeds. 

The objective is to derive analogous expressions of d{uj) 
and G for anisotropic elastic materials. This will be 
achieved by explicit calculation of the admittance ten- 
sor A, defined in Section |TT1 combined with a general 
relation between d(tt'), G and A. The spectral density 
and the participation tensor in the presence of material 
anisotropy do not appear to have received much atten- 
tion. Some work on the related issue of admittance in 
bounded anisotropic thin plate systems has appeared^. 
Weaver— considered isotropic plates of finite thickness 
and infinite lateral extent. Tewary et ali^ derived an 
expression for the admittance at the free surface of an 
anisotropic half space as a double integral. Here the fo- 
cus is on infinite systems, and the modal density per unit 
volume in this limit. Finite structures, such as plates 
both thin and of flnitc thickness, will be considered in a 
separate paper. 

Our principal results arc that the modal spectral den- 
sity per unit volume and the participation tensor are 
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TABLE I. The form of the participation tensor G for the dif- 
ferent material symmetries. TI, tet and trig are abbreviations 
for transverse isotropy, tetragonal and trigonal symmetries, 
respectively. The e unit vectors are defined by the symmetry, 
while a, b and c result from averaging. The positive numbers 
a, P and 7 are constrained as indicated in order to satisfy Eq. 
©■ 



Material symmetry 
isotropic, cubic 
TI, tet, trig a + 2/3 = 3 
orthotropic a + /3 + 7 = 3 
monoclinic a + /3 + 7 = 3 
triclinic a + /3 + 7 = 3 



I 

Qe (g) e + /3(I — e (g) e) 
aei (g) ei + /3e2 ® 62 + 763 ® 63 
ae (g) e + /3a ® a + 7b ® b 
aa ® a + /3b ® b + 7c ® c 



given by 



d{uj) 



(trQ" 



G =3- 



-3/2\ 



(trQ" 



(4a) 
(4b) 



where Q(n) is the acoustical or Christoffel tensor for 
plane waves propagating in the direction n, and (/) is 
the orientation average of a function that depends on the 
direction, 



47r 



(5) 



In an isotropic solid (|4a|) reduces to ([3]) and G is sim- 
ply the identity I. After deriving the remainder of 
the paper will explore its implications, in particular the 
form of G is investigated, and the parameters in Table 
I deduced. It is interesting to note that the material 
constant that determines the density of states of diffuse 
waves, tr(Q~'^/^), also defines the Debye temperature 6 
of a crystal. Thus (see Chapter 9 of ReffTT 



k\VatT{Q 



-3/2\ 



1/3 



(6) 



where h is Planck's constant, k is Boltzmann's constant, 
and Va is the volume per atom or lattice site. Fedoro\»i^ 
provides a detailed discussion of tr(Q~'^/^) in this con- 
text. The emphasis in this paper is on the more gen- 
eral tensor (Q '^^^) although connections with Fedorov's 
analysis will be mentioned later. 

The outline of the paper is as follows. The admittance 
tensor A is defined and calculated in Section [Til from 
which the main result (|4|) follows. Several alternative 
representations of the fundamental quantity Q^^^^ are 



developed in Section IIIII In particular it is shown that 
G for transverse isotropy can be evaluated as a single 
integral. Weak anisotropy is considered in Section IIVI 
and numerical examples are presented in Section fVl 



II. DERIVATION OF d AND G 

A. Admittance tensor 

The admittance A is a second order tensor defined by 
the average power radiated by a time harmonic point 
force F according to 



n F • A • F. 



(7) 



Alternatively, A is equal to the power expended at the 
source point - which is the more conventional definition of 
admittance, as the the inverse of drive point impedance. 
The admittance is clearly related to the auto-correlation 
of the Green's function, and as such is a special case of 
the two-point cross correlation of the Green's function^. 
The important connection for the present purposes is the 
relation between the radiation from a point force and 
the diffuse wave density^i^. In the present notation this 
becomes 



A = d(uj)G. 

Up ^ ' 



(8) 



A short derivation of is given in Appendix The 
admittance of isotropic bodies is simply determined from 
Eq. ([3]) and G = I. Our objective here is to calculate 
A for anisotropic solids, and then to use the result to 
determine d{uj) and G. 

The central result for A is the following: The second 
order symmetric admittance tensor of Eq. (O that deter- 
mines the total power radiated to infinity from the point 
source averaged over a period, is 



where Q(n) is the acoustical tensor. 



Qik{n) = djkirijni with Cijki = -CkUj- 



(9) 



(10) 



The elastic moduli (stiffness) Cijki have the symmetries 



C. 



jikl : 



and thus have at most 



Cijki — Ckiij and Cijki 
21 independent elements. Note that A has dimensions 
of admittance (inverse impedance). We next derive Eq. 
^ by explicitly calculating the admittance for a time 
harmonic point force. 



B. Radiation from a point force 

The displacement resulting from a point force F cos ujt 
at the origin is u(x, t) = Re u(x, cj)e~*"* where u satis- 
fies 

Cijkiuk.ji + puj'^Ui = -Fi(5(x), -00 < xi,X2,X3 < 00. 

Here p is the mass density and (5(x) is the three- 
dimensional Dirac delta function. The equation of mo- 
tion may be written 



Q(V)u-f tj2„ ^ _i^(x)F 



(11) 
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and the problem definition is completed by the require- 
ment that the energy radiates away from the point 
source. 

The solution to pT|) in a solid of infinite extent is well 
known. For our purpose we will find the following repre- 
sentation from Norrisi^ (Eq. (3.22)) useful for determin- 
ing the admittance: 




Here Ai,A2,A3 are the eigenvalues and qj,q2,q3 the 
eigenvectors of Q(n), which then has the spectral de- 
composition 

Q(n) = AiQi ® -I- A2q2 + Agqg (g) qg. (13) 

1 /2 

Also, kj = w/Aj are the wavenumbers of the three 
distinct branches of the slowness surface defined by the 
eigenvectors. The first integral in (jl2p is around the unit 
circle formed by the intersection of the plane n • x = 
with the unit n— sphere. This is just the static Green's 
function of elasticityiS. The important dynamic quan- 
tity is the second integral which is evaluated over the 
sphere {|n| — 1}. In order to make this more apparent, 
we rewrite ([T^ as 

u^u\..o + ^t{e-'-^)F, (14) 

and note for future reference that the first term on the 
right hand side is real valued. 

The average power radiated per period is equal to the 
power expended by the force 

n=lim— / dt coswtF • v(0,t), (15) 
x^o 27r 7 ^ ' ^ ^ 



where v(x, i) = Re ( — iwu(x, a;)e~"^') is the particle 
velocity. Thus, 

The spectral decomposition ([T3| implies that 

which together with Eq. ^ proves the main result 

The scalar d{u!) and the tensor G are defined such 
that their product is 12p/7r times the admittance A, see 
Eqs. dSl), (HI), dH) and ^. This defines d and G to 
within a constant, which is determined uniquely by the 
constraint tr G = 3. We therefore obtain the general 
results of Eq. As discussed, d is the generalization 
of the classical density of states per unit volume, ^ for 



isotropic solids, and the participation factor tensor G 
describes the directional distribution of the energy at a 
point. While it is convenient to consider them separately, 
d and G are both defined by the averaged tensor (Q~^^^) , 
which will be the focus of the remainder of the paper. 

Before considering the properties of d and G we note 
that the isotropic modal density of states follows imme- 
diately from (|^a|l . Starting with the acoustical tensor for 
an isotropic solid, 

Q(n) = c^n (X) n + (I — n (K) n), isotropy, (17) 

we have Q"'^/^ = cj'^n® n + c^^{I — n® n). Then using 
the fact that (n (g) n) = it follows that 

(Q-3/2) = i(cr^ + 2cr^)i. (18) 

Hence, the density of states per unit volume is d = 
^^(c;""^ -|- 2c^^)~^, in agreement with the well known 
identity and G = I, as expected. 

III. Q^/^ p^f^Q RELATED QUANTITIES 

The key quantity is the tensor Q^'^''^ and its direc- 
tional average. In practice, this may be evaluated nu- 
merically without difficulty. It is however useful to ex- 
amine semi-explicit forms for the tensor, both for general 
anisotropy and for specific symmetries, particularly the 
case of transverse isotropy. We begin with two alternative 
and general formulations based on the spectral properties 
and the invariants of the acoustical tensor. 



A. General representations for arbitrary anisotropy 

1. A method based on invariants 

Functions of a positive definite tensor can be simplified 
using the Cay ley-Hamilton formula for the tensor, which 
for Q is 

Q^-/iQ' + /2Q-/3l = 0. (19) 
The principal positive invariants of Q are 

/i=trQ, /2 = i(trQ)2-itrQ2, /g^detQ. (20) 

Based on these fundamental properties, it can be shown 
that 

=[(^1^1^ + ,^^3/2 + z2/3)(Q2 - /iQ + /2I) 

+ iihhiQ ' hi) - l!i\/[{iii2 - (21) 

1/2 

where ii, «2 and 13 arc the positive invariants of Q ' 
which can be expressed as functions of the invariants Ii 
and /g, see below. Details of the derivation of (pij) are 
given in Appendix B. 

The appeahng feature of Eq. for Q"^/^(n) is that 
it only involves powers of Q, its three invariants, and the 
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additional invariants 11,^2 and 13. These are related to 
h, h and/3 hy^ 

il-2i2^h, il-2iii3^l2, il^h- (22) 
1/2 

The last implies = , while expressions for ii and 12 
are given by Hoger and Carlsonii and by Norrisi^. For 
instancei^, 



t2 =\ll2-hlP + 2./hP + Vhl^. 
where /3 is any eigenvalue of Q, e.g. 

/3 + [(e + ^e-{i!-3i2r] 



(23a) 

(23b) 
(23c) 



1/3 



+ [i^-^je-ii!-3i2r] 

^=i(2/3- 9/1/2 + 27/3). 

Note thatH i^i^ - 43 = dct(iil - Q^^^) > 0. 
Taking the trace of Eq. (j2ip gives 



(24a) 
(24b) 



, P,-3/2 _ ih + ^2)/2/3 + (/| - 2/l/3)n^3 ~ 3/| 
(«1«2 - «3)/| 

(25) 

This quantity, when averaged over all orientations, gives 
the density of states function d{ijj) of Eq. ((4a|) . Hence d 
can be calculated from the invariants Q and the derived 
invariants zi, 12, is- 



2. A spectral representation 

The second form for Q^'^''^ is based on the spectral 
decomposition pT]) . The latter can be expressed in a 
form that docs not explicitly involve the eigenvectors. 



Q"'/' = A7'/'N(Ai) + Ar''N(A2) + A3"'''^N(A3). (26) 

The second order tensors N(Aj), which are alternative 
expressions for the dyadics formed by the eigenvectors, 
% ^ I can be expressed in terms of Q using Sylvester's 
formula 



,-3/2 



-3/2, 



N(A,n) 



AQ^ + (A-/i)AQ + /3l 

A3 + (A-/i)A2 + /3 ■ 



(27) 



The identity (|26p is derived in Appendix IB] 

Calculation of (pS)) requires knowledge of the three 
eigenvalues, which arc zeros of the characteristic poly- 
nomial defined by Eq. p9)) . 



(28) 



p(A) = A-^-/iA2 + /2A-/3 



The eigenvalues {Ai, A2, A3} can be expressed in terms of 
the invariants as 



1 



1 



{/3, -(/i - /?) ± - V(/i - /3)2 - 4/3//?}, (29) 



where f3 is defined in (|24a[) . Everyii derived alternate 
closed-form expressions based on the trigonometric solu- 
tion of the characteristic cubic. The alternative version 
of Eq. 



(30) 



which is the starting point for Fcdorov's calculationii of 
the trace. 



B. Transverse isotropy 

Transverse isotropy or hexagonal symmetry is an im- 
portant class of anisotropy. It occurs in many practical 
circumstances, whether from layering in the earth to lam- 
inated composite materials, or from underlying crystal 
structure. It is the highest symmetry for which the par- 
ticipation factor tensor is not the identity, since G = I 
under isotropy and cubic material symmetry. Wc now 
demonstrate that the evaluation of d and G may be re- 
duced to the evaluation of two single integrals, one for 
(trQ~'^/^) and one for the parameter a that defines G, 
see Table I. 

Transversely isotropic solids have five independent 
moduh: cn = C22, C33, C12, C13 = C23, C44 = C55, 
cgg = ^(cii — C12). Let e be the axis of symmetry. The 
SH slowness decouples to give 

Q = A3(n-e)q3®q3-|-Q_L, (31) 

wherein (p. 95) 

A3(n • e) = C66 + (c44 - C66)(n • e)^, (32) 

and q3 = eAn/|eAn|. The 2-dimcnsional symmetric 
tensor isf^ 

=[c44 + (C33 - C44)(n • e)^]e (g) e 
+ [cii + (c44 - cii)(n • e)2]d (g d 



-I- (ci3 -I- C44)n • e\/l - (n • e)2[d (g) e + e (g) d] , 

where d = e A q3. Replacing n • e by the integration 
parameter ^, it follows that 



(A; 



-3/2 



1 



q3 ® ^3) 



(33) 



where Ij^ projects onto the plane perpendicular to e, 

I±=I-e®e. (34) 

3 /2 

It remains to consider the orientational average of Qj^ 

The tensor Qj^ satisfies a quadratic Cayley-Hamilton 
equation 



Qi- JiQ± + J2l± = 0, 



(35) 



with Ji = trQj^ = Ai -I- A2 and J2 = det Qj_ = A1A2. 
Similarly, the Cayley-Hamilton equation for the square 
root is 



(Qf)^-jiQf 



.72 1 J 



0, 



(36) 
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where ji = tr and j2 = det satisfy Ji = jf — 
2j2, J2 = il' ^i^d are therefore related to Ji and J2 by 

ji = y^Ji + 2VJ2, j2 = V^. Using Eqs. ([33 and ([Ml), 
respectively, leads to the identities 

Ql' =J2'' M - J2)I± - JiQa] , (37a) 
Qy'=jr'(Q±+J2l±). (37b) 
Multiplication of these and further use of ([35]) leads to 

q;3/2 = ^L_[(j^ _ Q ) _ (38) 

Again using ^ = n • e, we have 

(trQ-3/2^ = 
1 

J de [j2"'/'(Ji -^2)\/ji + 2v/j2 + A-'/'io] , 





(trQ-3/2) 
a 



where 

a = Cii+C44, fo = C33-Cii, 

C=C44-Cii, d=CiiC44, 

e = C11C33 - Ci3 - 2c44(cii + C13), 

/ = -C11C33 + Ci3 + C44(cii + C33 + 2ci3). 

IV. WEAK ANISOTROPY 

Although the general expressions for the modal density 
d and the participation tensor G are not difficult to com- 
pute, it is often the case that the medium is to a first ap- 
proximation isotropic, and appropriate approximations 
can be made. The state of small or weak anisotropy is 
defined relative to a background isotropic medium, and it 
is important to select the latter properly. In this Section 
we calculate d and G in the presence of weak anisotropy. 
Fedorovii provides a detailed analysis of the expansion 
of tr(Q~^^^) to arbitrary orders in the perturbation pa- 
rameter. Our emphasis is more on obtaining estimates of 
the tensor (Q"'^/^), which is not discussed explicitly by 
Fcdorov. Wc begin with a description of the comparison 
isotropic moduli and then proceed to calculate the first 
two terms in a perturbation series for d and G. 



and from Table I, 



3 } {Ji + v7^){Ji - e ■ ■ e) - J2 



The modal density parameter (trQ""^/^) and the scalar 
a that defines the participation tensor can therefore be 
expressed as single integrals, which follow from the above 
results and Eqs. ([?T|) through as 



A. Background isotropic moduli 

Regardless of the level of the anisotropy it is al- 
ways possible to define a unique set of isotropic moduli 
which minimize the Euclidean distance between the ex- 
act set of moduli and the equivalent isotropic moduli^^. 
This procedure is equivalent to requiring that the mean 
square Euclidean difference in the slowness surfaces is 
minimali^iii. Thus, let the background isotropic moduli 
be 

4%i ^ (^f^ij^ki + CfiS^kSji + SiiSjk - 2dijSki), (40) 

where ci and ct are the effective longitudinal and trans- 
verse wave speeds. These are defined by simultaneously 
minimizing the quantity (|Q — QqP) with respect to both 
ci and Ct, where Qo(n) is defined by the moduh c^^li- The 
unique solution is 

= ^ tr Ci, c? = i tr Ct, (41) 
where the second order tensors of reduced moduli are 
2 1 r ~ 1. -1. 

(42) 

The background Lame moduli A and fi are obtained using 
cf = (A + 2/i)/p and Cj = /i/p. The elements of C; and 



{a + be - Vd + ee+fe)Ja + be + 2^d + ee + fC 



{d + ee + le?'^ 



[c66 + (C44-C66)C']3/2_ 



(trQ" 



-3/2\ 



(cn + ce){a + be + ^d+ee + fe) -{d+ ee + fe) 



(d + ee + feY'^\ja + be + 2^d + ee + fe 

I 



(39a) 
(39b) 
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Ct follow from 

^Cii + Ci2 + Ci3 Ci6 + C26 + C36 C15 + C25 + C35 ' 

CjjfeA: = I C16 + C26 + C36 C12 + C22 + C23 C14 + C24 + C34 

yCi5 + C25 + C35 Ci4 + C24 + C34 Ci3 + C23 + C33^ 

^Cii + C55 + C66 C16 + C26 + C45 Ci5 + C46 + 035^ 

Cifcjfe = I C16 + C26 + C45 C22 + C44 + C66 C24 + C34 + C56 

, Ci5 + C46 + C35 C24 + C34 + C56 C33 + C44 + C55 , 



B. Perturbation analysis 

Let 



(43) 



where the nondimensional parameter e is introduced only 
to simplify the perturbation analysis. In practice e is set 
to unity on the assumption that the additional moduli 
'^^^^ are small in comparison with the isotropic 



Ciikl 



background. 

We seek expansions in powers of the small parameter 
e. The key quantity Q~'^^^ will be determined as the 
product of Q^^ and Q^^^. Based on (^51) . the acoustical 
tensor is 

Q = Qo + £Qi> (44) 
and simple perturbation gives 

Q"' = Q(7' - e(Qf7'QiQo"' + Qo^'QiQo ') + 0{e^). 

Let 

then Si satisfies 

(45) 



Qo^ Si + SiQq^ — Qi- 



In order to calculate Q and also the square root of 
Q, we now use the fact that the leading order moduli c|^],j 

1 /2 

are isotropic. The explicit form of Qg' follows from Eq. 
P?)) and the identity 

= c2™n ® n + c?'"?, (46) 

where m is any real number and P = I — n (g) n. Equa- 
tion can be solved by observing that Qi may be 

partitioned = Q*^^ + Q^' + ' where Q^^^ = 

n • Qi • n)n ® n, ^ = PQiP and Qf ' = PQi • 

nig)n-|-n(g)PQi - n. Assuming a solution of the form 

Si = PiQi^^ +P2Qi'^^ + P^Q^i \ the coefficients can be 



^ Qf^ (47) 



determined easily from Eq. p5|) . i.c 

Si = ^Qi^' + ^Qf ^ + -- 

2ci 2ct ci + Ct 



Combining the asymptotic expansions for Q ^ and Q^^^ 
gives 

Q-3/2 = Q-3/2+£Vi+0(£2), (48) 



where 

Vi =Qo ^Si - Qq ^QiQq — Qg ^QiQo 
-„An r (c? + c? + QCt) 3. 

X [Qi • n (K) n + n (K) Qi • n] 

r(cf +cf + QCf) 3 3 T , ^ 

The orientational average can then be effected 

using the identities 



ijkl 1 



(n^njnkninpnq) ^-{StjKkipq + 6tkKjipg+ 



S,iK, 



il-^kjpq ^ip-^^ kljq ^~ ^iq-^klpj 

The resulting expressions for (Q^"^^^) is 



3 V? 



2^5 ikjk 



2 . {cf + cj+QCt) 3 . (1) (1) 
15 L cfc!ici+ct) 2c5 J ^ "^ ^'^-■'"'^^ 

^ r2 (£L±fL±£i£t) _ J ^1 

105 L cfcUci + Ct) 2cf 2cf5J 

X [^.,(ci:;i, + 2ci:;L) + 4(c«, + 2cW,)] j + 0{e') 



We note that both ci^], and c,-^] ■ vanish by virtue of 
the choice of the background isotropic moduli. This im- 
plies that the trace of (Q '^''^) differs from the isotropic 
approximant only at the second order of anisotropic per- 
turbation, 



tr(Q- 



-3/2\ 



2 1 , o 

:j + 3 + o(^ 



(49) 



This is in agreement with Fedoro\^ who also provides 
explicit forms for the higher order terms; for instance, 
the expansion for cubic crystals up to fourth order in 
the perturbation is given by Eqs. (50.12) - (50.14) of 
Refjrg. The leading order approximation of Eq. 
when combined with the identity (j4b[) . gives 



/ 2 1 \-i J 3 (1) 

i_ ro (c' + Ct +QCt) , 1 _ 1] 
35 L cfcf{ci+ct) cf c^J 

x(cyL+2c«,)|+0(e^). 

Ignoring terms of order and then setting £ — > 1 yields 
the leading order approximation to the participation ten- 
sor as 

G « I + ai{I - c^^Ci) + at{I - q-'Q), (50) 
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FIG. 1. The non-dimensional parameters a; and at as a func- 
tion of the Poisson's ratio v. 



FIG. 2. The non-dimensional parameters ai, 02 and 03 for 
weak transverse isotropy as a function of the Poisson's ratio 
V of the background medium. 



where the non-dimensional coefficients are 



ai 
at 
and 



7(2-f K-3) V3 
3 



Ct " 



(51a) 
(51b) 

(52) 



Figure [T] shows ai and at as functions of the Poisson's ra- 
tio V, using = 2(1 - i/)/(l - 2v). Note that 1.27. . . < 
at < 3/2 for < < 1/2 while a/ « -^(1 - 2iy)~^ as 
ly 1/2. 



C. Transversely isotropic materials 

As an example of the general perturbation approach, 
we consider the particular case of TI materials. We take 
the axis of symmetry (e in Section|TTT| in the 3— direction, 
so that 



'C11+C12 + C13 

Cii -I- C12 -I- Ci3 

C33 + 2ci3^ 

' Cii + C44 -I- C66 

Cii + C44 + C66 

C33-f2c44, 



where cge = 5(^11 ~ C12). The wave speeds in the back- 
groimd isotropic medium are then, 

= ^(8cii +3c33 + 4ci3 + 8c44), (53a) 
= ^(2cii + 2c33 - 4ci3 + 12c44 + lOcee). (53b) 



According to Table I the participation tensor is defined 
by a single parameter, a, which to leading order is unity. 
Let 



so that 



G = 



a = 1-2/3, 



'1 + /3 
1 + /3 
l-2/3y 



(54) 



(55) 



Applying the general perturbation theory we find that 
the leading order correction to the isotropic participation 
tensor is given by 



ai 



(-4cii -I- 3C33 -I- Ci3 -I- 2C44) 



at 

met 



(-C11 -I- 2C33 - Ci3 -I- 3C44 - 5C66), (56) 



where ai and at are defined in l|51ap . 

Thomson's anisotropy parameters'^ e, 7, S provide a 
means to characterize weakly anisotropic TI materi- 
als. The parameters are defined e = (cn — C33)/(2c33), 



S = [(Ci3 + C44) 



(C33 



C44) ]/[2c33(c33 - C44)], 7 



(c66 — C44)/(2c44), and are commonly used in geophysical 
applications to describe rock properties. The correction 
term /3 can be expressed in terms of the Thomsen param- 
eters as, 



/? « aie + a2S + a^j, 
where the coefficients ai, 02 and 03 are 



ai 



15 



a2 



ai K at 



15 30 ' 



aa 



at 



(57) 



(58) 
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TABLE II. The participation matrix G for a variety of 
anisotropic materials. Sym denotes material symmetry: 
transversely isotropic (TI), tetragonal (Tet) or orthotropic 
(Orth). The Frobenius (p=2) norm is used to compare G 
with the isotropic result (I) and with the perturbation ap- 
proximation G defined by Eq. (|50|) . dist is a non-dimensional 
and invariant measure of the anisotropjsi^, equal to zero for 
isotropy. dist> 1 signifies considerable anisotropy. 



Material 


Sym 


Gil 


G22 


G33 


|G-I| 


|G-G| 


dist 


Beryllium" 


TI 


1.05 


1.05 


0.89 


0.13 


0.00 


0.22 


Sulphur" 


Ort 


0.95 


1.32 


0.73 


0.42 


0.11 


0.95 


Cadmium" 


TI 


0.73 


0.73 


1.55 


0.67 


0.10 


1.02 


Barium titanate'' 


Tet 


0.81 


0.81 


1.39 


0.48 


0.01 


1.11 


Rochelle salt" 


Ort 


1.38 


0.65 


0.97 


0.52 


0.09 


1.16 


Zinc" 


TI 


0.71 


0.71 


1.58 


0.71 


0.14 


1.17 


Graphite/Epoxy'^ 


TI 


1.38 


1.38 


0.25 


0.92 


0.81 


2.35 


Tellurium dioxide'' 


Tet 


1.30 


1.30 


0.40 


0.74 


0.72 


2.87 


Mercurous iodide"* 


Tet 


1.37 


1.37 


0.26 


0.91 


0.14 


3.02 


Spruce" 


Ort 


1.35 


1.63 


0.02 


1.22 


1.30 


5.59 



"Elastic moduli from RefHsl. 
''From Refill 
Trom Refill 
"^From ReflU 



V. EXAMPLES AND DISCUSSION 

The participation matrix was computed for many 
anisotropic solids. Table II summarizes the results for a 
selection of materials with anisotropy ranging from weak 
to strong. The table provides the numerical values of 
diagonal elements of G (there are no off-diagonal ele- 
ments for the symmetries considered). In each case the 
elements sum to three, Gn + G22 + G33 = 3, although 
the individual numbers can differ markedly from unity. 

In order to quantify the level of anisotropy, the table 
also shows the number dist. This is a nondimensional 
positive measure of the degree of anisotropy of a set of 
anisotropic elastic constants, dist is chosen here as the 
log-Euclidean distance or length from isotrop y al- 
though other measures are possible, see Norris^^ for a 
comparative discussion. The log-Euclidean distance has 
the advantage that it is invariant regardless of whether 
the compliance or stiffness tensor are considered. We use 
dist as a convenient and simple measure of the degree 
of anisotropy. Appendix [C] provides a little more detail 
on its exact definition, including a short Matlab script to 
compute dist. 

Large deviations from the isotropic participation ten- 
sor are apparent. Consider the ratio R of the largest 
to smallest element of G. Even for small to mod- 
erate anisotropy, such as Cadmium we see that R = 
G33/G11 > 2. The ratio becomes much larger for the 
more anisotropic materials considered. Spruce is in- 
cluded because of its enormous ratio, R 80. These 
ratios can be compared with the results for the rela- 
tive partition of the diffuse wave energy at the free sur- 
face of an isotropic solid. If 63 is the normal to the 
surface, then the calculations of Weaver— indicate that 



1 < G33/G11 < 1.25 where the lower (upper) bound 
is reached as i> approaches 1/2 (0). The upper bound 
« 1.25 is approximate and based on Fig. 3 of RefH. 

The numbers in Table II indicate that the perturba- 
tion approximation is adequate for small anisotropy. This 
can be characterized loosely as < dist< 1, and strong 
anisotropy is dist> 2, roughly. The examples in the Table 
suggest that the weak anisotropy approximation is not 
useful in the presence of strong anisotropy. This is evi- 
dent from the fact that the errors |G— 1| and |G— G| are of 
the same order of magnitude for the strongly anisotropic 
materials, whereas |G — G| is much less than |G — 1| for 
weak anisotropy. 

We note that for all materials considered the numerical 
calculations show Eq. underestimating tr(Q~^^^). 

However, the more refined perturbation expansion of 
tr(Q~^^^) by Fedoro\»ii suggests that this is not a uni- 
versal result. 

The dependence of G and d(uj) on the moduli is ob- 
viously complicated by virtue of the averages required 
in Eq. ([4]). However, the formula ((50)l for G for weak 
anisotropy illustrates the dependence more explicitly. 
The form of the matrices C; and Ct imply that only 12 
combinations of the 21 independent anisotropic moduli 
enter into the first term in the perturbation expansion. 
For orthotropic materials, with 9 independent moduli, 
this number reduces to 6, and the matrices C; and Ct 
arc then diagonal. In the case of weak TI only two com- 
binations of moduli infiuence G, see Eq. ([56|. 

The non-dimensional tensor G also has important im- 
plications for radiation from a point source. The connec- 
tion follows from the relation ([8]) between G and A, com- 
bined with the correspondence between the drive point 
admittance tensor and the radiation efficiency in Eq. ([7|. 
Thus, the direction in which a force must be applied to 
most efficiently radiate power is the principal direction of 
G with the largest element. Conversely, the least amount 
of power is radiated if the force is directed along the prin- 
cipal direction with the smallest element. For instance. 
Table II indicates that a point force of given magnitude 
will radiate most power in Cadmium if the force is di- 
rected along the axis of hexagonal symmetry. The situa- 
tion is reversed for aligned graphite/epoxy, where forcing 
along the fiber direction produces the least amount of to- 
tal radiated power. 

The inverse problem of determining anisotropy from 
measurements of G is clearly ill-posed. However, pos- 
sible measurement could be advantageous in particular 
circumstances. Consider for instance, 3-component mea- 
surement of the displacement downhole in a borehole 
environment. Assuming the frequency is such that the 
wavelengths are large compared with the bore radius, 
the 3-component data is sufficient to compute the auto- 
correlation and hence G. The principal directions of G 
and the relative magnitude of its diagonal elements pro- 
vides significant information about the local geostratig- 
raphy and formation properties. 
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VI. CONCLUSION 

We have derived general formulas for diffuse waves 
in anisotropic solids. The main results are concise ex- 
pressions for the modal density per unit volume and fre- 
quency, d(iu) of Eq. (|4ap . and the participation tensor G 
of Eq. (|4bp . The latter is a material constant with one 
or two independent constants, and with principal axes 
dictated by the material symmetry. In the absence of 
symmetry the participation tensor defines principal axes 
for diffuse wave energy distribution, and for radiation 
efficiency. Calculation of ditu) and G requires, in gen- 
eral, averaging over the surface of the unit sphere. Single 
integrals suffice for transverse isotropy, with the impor- 
tant quantities given in Eq. ([M]) . In the case of weak 
anisotropy, a perturbation scheme produces explicit for- 
mulas, Eqs. (Uni and (|5D|) . The main quantity in all 
cases is the second order avcrag ed tensor (Q^^/^). We 
have illustrated the results through calculations for sev- 
eral materials. These display the main effects that would 
occur in all anisotropic solids. In particular, the devia- 
tion G from the unit identity tensor can be significant. 
Ratios of 2 or more for the relative magnitude of diffuse 
wave energy in different directions in crystals can occur 
under moderate levels of anisotropy, with far larger ratios 
possible in realistic materials. 
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APPENDIX A: DERIVATION OF EQ. ^ 

We use an argument based on a modal representation^ 
for the solution to the point force problem, 

(p—^ — L)u ~ F6(x ~ xq) cosuit, (Al) 

where L is a second order differential operator. The re- 
sulting velocity v = du/dt may be found by standard 
means as 



1 V- 



-imF ■ Um(xo) Um(x) 

cjH, — — zO 



where the modes u,„(x)e are solutions of the ho- 

mogeneous equation jATJ, with the properties 

(5(x - xo)I = ^ u,„(x)u™(xo), 

rn 

/ dxu„i(x) • u„,(x) = 1. 
Jv 

The power output averaged over a cycle is therefore 

n(xo,w)= — / di coswi F • v(xo, i) 
Jo 



= -^[F.u,„(xo)]2Re 



(A2) 



The 



strict 



non-dissipative 



limit 



of 



Rc[— — — *0)^-^] is 7ra;(5(w^ — oj^) = 
^Tr6{ujm — uj) where S is the Dirac delta function. 
However, modal overlap in the presence of non-zero dissi- 
pation spreads the influence over many modes. The effect 
is to make Re[—iuj{uj^ — — iO)^^] ^TTf{uJm — i^) 
where f{v) is smooth with bounded support in 
€ {— fi, say, and unit sum: 



X! /(^™ -oj) = 1. 



(A3) 



Here indicates the sum over modal frequencies 

ujm G {w — ri,cij + Q}. Using the density of modes, 
Vd{u!m), to replace the sum over modes in (jA2p by a 
sum over modal frequencies, gives 

ttV , ' ^ 

H(xo,cj) = — ^d(w,„)/Kn-t^)[F-u,„(xo)]^ (A4) 

We now make the assumption that the support of /(i^) 
is small enough that the modal density function, d(ujjn): 
may be replaced by d{uj). This is perfectly reasonable 
based on known forms for d{uj), e.g. Eq. At the same 
time, we assume that the support of f{h') is sufficiently 
large that we may use the equipartition of energy among 
modes to make the replacement (see Eq. ([T])) 



y^/(w™ - w)u„ y^u® u = -G. (A5) 

Hence, 



H(xo,c^)-— d(co)F.G.F, (A6) 
12p 

and since F is arbitrary, the admittance A follows from 
the definition of H in ([7]). This completes the derivation 
of the identity ©. 



APPENDIX B: DERIVATION OF EQS. (gl]) AND ([26ll 

The Cayley-Hamilton relation for Q is p(Q) = 0, where 
p is the characteristic cubic polynomial defined in Eq. 
(PS)) . and /i(n), /2(n), ^3(n) are the invariants define in 
Eq. Thus, 



h — -^1 +-^2 + ^3, I2 — A1A2 + A2A3 + A3A1, /3 — A1A2A3, 

and since Aq = it follows that the invariants are all 
positive, /i > 0, /2 > and I3 > 0. Multiplying by 
Q""'^ and yields equations for the same quantities: 



3 



Eliminating Q^^ gives an equation for Q^^: 

Q-2 = 73-2 _ (7^/2 _ /3)Q + (/| _ 7^/3)1] 



(Bla) 
(Bib) 
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1/2 

We next derive a similar type of equation for Q ' using a 
method due to Hoger and Carlsonii. The product of this 
with Q^^, combined with the Caylcy-Hamilton equation 
(fT9)) yields the desired relation (PT|) . 
First we note the general expression 

(Q-AI)-i = 

^ [ - Q2 + (/i - A)Q - (A2 - /lA + h)!] , (B2) 
p{X) 

where p is the characteristic polynomial for Q, from Eq. 
(PS)) . The identity (jB2[) may be checked by direct mul- 
tiplication and use of Eq. p9p . The square root tensor 
R = Q ' satisfies the Caylcy-Hamilton equation 

R^-iiR^+i2R-i3l = 0, (B3) 

where ii, 12 and 13 arc related to the invariants of Q by 

h=il- 2i2, h=il- 2^l^3, h = «3- (B4) 



Explicit formulae for ii, 12 and 13 are given in (|23ap . 
Rearranging (|B3|) as R(R^ -1-12!) — «iR^ -I- 13! and using 



R" 



Q gives 



R= (nQ + «3l)(Q- 



«2J 



(B5) 



Application of (|B2p along with some simplifications using 
(jB4p . such as p{—i2) = —(is — *i*2)^! yields 



Combining Eqs. ((B2l) and ((B6)l gives Eq. ([211) • Alterna- 
tively, 



where 



Q 



-3/2 



aQ^ + bQ + cl, 



(B7) 



6 = 



/|(i3-ilZ2) 

hh{i\ - 12) + (/1/2 - /3)n»3 

/|(i3-Hi2) 
/| + 12/3(^2 - ij) + (/1/3 - iDnh 

^3(«3-«l«2) 



(B8) 



The second form (|26|) for Q "^^^ is based on the identity 
(fT7)) . The tensor products of eigenvectors for A^ satisfy 



(Q-A,I)(Q-AfcI) 

(A,- A,)(A,-Afe) 



i ^ i ^ k ^ i (no sum). 



This follows, for example, by eliminating the other two 
tensor products using the spectral expressions for I, Q 
and Q^. The dependence on \j and A^ can be removed 
in favor of Aj and the invariants Ii and I^, and hence Eq. 
([27]) . Note that the latter can be expressed 

N(A, n) = [AQ2 + (A - /i ) AQ + /3I] , (B9) 

where p'{x) is the derivative of the characteristic poly- 
nomial. This indicates that the general expression ([77]) 
is invalid at double roots where the slowness surface ex- 
hibits degeneracy, and proper limits arc required. The 
possibility of such points does not present a practical im- 
pediment to numerical integration. 



APPENDIX C: THE LOG-EUCLIDEAN DISTANCE 

The procedure^ is to first calculate an effective 
isotropic set of moduli analogous to c[°],; of Eq. ([^0]) 
but for the matrix logarithm of the 6-dimensional Voigt 
matrix of moduli C/j. Some matrix factors are required 
to convert from the Voigt notation. The following Matlab 
lines compute dist if C is the 6x6 Voigt matrix. 

J = 1/3* [1 1 1 0] 1 1 0] ; 
K = eye(6)-J; 

T = diag([ 1 1 1 sqrt(2)*[l 1 1] ]); 
L = logm(T*C*T) ; 

dist = norm(logm( J*exp(trace(J*L)) 

+ K*exp(l/5* trace(K*L)) )- L ,'fro'); 
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